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We report density functional theory calculations for Mo3S7(dmit)3. We derive an ab initio tight- 
binding model from overlaps of Wannier orbitals; finding a layered model with interlayer hopping 
terms ~ 3/4 the size of the in-plane terms. The in-plane Hamiltonian interpolates the kagome 
and honeycomb lattices. It supports states localized to dodecahedral rings within the plane, which 
populate one-dimensional (ID) bands and lead to a quasi-ID spin-one model on a layered honeycomb 
lattice once interactions are included. Two lines of Dirac cones also cross the Fermi energy. 


Understanding emergence, the observation that sys¬ 
tems can have fundamentally different properties from 
their constituent parts, is one of the most important 
challenges in modern physics [1-6]. For example, spin 
liquids may realise emergent gauge fields and deconfined 
fractionalised excitations [7]. However, predicting emer¬ 
gent behaviour in condensed matter systems is incredibly 
challenging. Thus, identifying motiefs that lead to cer¬ 
tain emergent properties is an extremely valuable step 
towards the goal of rational design of complex materials. 

Molecular crystals are an exciting test ground for find¬ 
ing and understanding new emergent states of matter. 
They often display competition between multilple emer¬ 
gent strongly correlated ground states [8, 9]. This, and 
the flexibility of organic chemistry, means that the emer¬ 
gent physics is often tunable by subtle chemical and phys¬ 
ical modification [9-11]. 

Mo 3 S 7 (dmit )3 is a single component organometallic 
molecular crystal that was first studied as a potential or¬ 
ganic metal. Experimentally, one sees an insulator with 
well-defined localised magnetic moments, but no long- 
range magnetic order down to the lowest temperatures 
measured (T = 2 K, compared to the magnetic interac¬ 
tion J ~ 100 K) [12]. Hence this material may have an 
emergent spin-liquid ground state [13, 14], just one of a 
few materials have been identified as candidate spin liq¬ 
uids [7, 8, 15, 16]. Fig. 1(a) shows the triangular shape 
of the molecule and the packing of molecules into a cor¬ 
rugated honeycomb lattice with a P3 (C 3l ) space group. 
The two molecules per unit cell are displaced along the 
c-axis by only slightly less than half the unit cell height, 
implying a subtle layering of the structure. 

In this manuscript, we report ab initio electronic struc¬ 
ture calculations on Mo 3 S 7 (dmit) 3 , which predicts lines 
of Dirac fermions and the emergence of ID behaviour 
from 3D electronic and crystal structures. We construct 
tight-binding models by producing a localized Wannier 
orbital (WO) representation of the low-energy bands. We 
then focus on the properties of the simplest such model, 
which allows us to explain the key features of the elec¬ 
tronic struture. 

In previous density functional theory (DFT) calcula¬ 


tions [12] antiferromagnetic order was found to give rise 
to an insulating state. The band structure predicted for 
the antiferromagnetic DFT solution is strongly ID, which 
led to proposals of quasi-ID models for Mo 3 S 7 (dmit) 3 
[12, 13], and proposals for spin-liquid states on the basis 
of these models [13, 14]. However, such long range an- 
tiferromagnetism has not been observed experimentally; 
so here we study the paramagnetic solution. Our param- 
eterisation therefore provides a suitable input for quan¬ 
tum many-body calculations that capture correlation ef¬ 
fects beyond those included in DFT. Here we construct 
a tight-binding model from first principles; we find that 
the model is fundamentally 3D as the largest in-plane and 
out-of-plane hopping integrals differ by less than a third. 
However, we show that interference effects produce two 
bands that are almost completely flat in-plane, while in 
the other bands the in-plane bands mass vanishes at the 
K and K’ high-symmetry points. We argue that a Hub¬ 
bard model at 2/3-filling on this lattice gives rise to a 
quasi-ID spin-one model. 

We show that both the flat bands (infinite mass 
fermions) and massless Dirac fermions occur in a 2D 
tight-binding model on a honeycomb lattice with each 
vertex decorated with triangles Fig. 1(b). We refer 
to this lattice as the ‘kagomene’ lattice as it smoothly 
interpolates between the kagome and graphene (honey¬ 
comb) lattices. Introducing an interlayer hopping be¬ 
tween stacked kagomene layers produces the ‘kagomite’ 
(kagome-graphite) lattice, which displays all the main 
qualitative features of our DFT band structure. Retain¬ 
ing only the three largest hopping integrals in the Wan¬ 
nier parameterization of the DFT band structure also 
produces the kagomite lattice. 

It has recently been shown that the exact ground state 
of the Kitaev model on the kagomene lattice is a chiral 
spin liquid [17]. The Kitaev model is a spin-1/2 model, 
corresponding to n = 3 electrons per triangular molecule 
[n = 4 in Mo 3 S 7 (dmit) 3 ]. It has been shown that, when 
spin-orbit coupling is included, the kagomene lattice sup¬ 
ports topological insulating phases and a quantum spin- 
Hall effect [18, 19]. Topological states have also been 
predicted on this lattice in the presence of staggered mag- 
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FIG. 1 : (Color online) (a) Crystal structure of Mo3S7(dmit)3 
projected along the c axis and (b) the kagomene lattice. In 
(a), one of the Wannier orbitals (WOs) for Mo3S7(dmit)3 is 
also shown. Note the predominantly dmit 7r-orbital character, 
with notable M03 d-orbital character. The other six WOs are 
related to this one by the C3 symmetry of the molecules and 
the inversion symmetry relating the two molecules in each 
unit cell. Sulfur atoms are shown in yellow, molybdenum in 
purple, and carbon in brown. In (b) two plaquette states, j 4 r 
and £>r/ , and a topologically non-trivial loop state are shown. 


netic fluxes [20]. 

We performed DFT calculations of the electronic struc¬ 
ture of Mo 3 S 7 (dmit )3 in an all-electron full-potential lo¬ 
cal orbital basis using the FPLO package [21]. The den¬ 
sity was converged ona (8x8x8) fc mesh using the PBE 
generalized gradient approximation [22]. Localized WOs 
were constructed from the six bands closest to the Fermi 
energy, and with these we computed real-space overlaps 



FIG. 2 : (Color online.) DFT and tight-binding band struc¬ 
tures for Mo3S7(dmit)3. The main panel shows the DFT 
band structure (solid black) along with bands from the tight- 
binding models constructed from the WOs with a large (r c = 
35 A, dashed green) and medium (r c = 20 A, dashed blue) 
real space hopping cutoffs. The r c = 35 A bands lie nearly ex¬ 
actly on top of the DFT bands. Lower inset shows the bands 
in a larger energy window, highlighting the gaps between the 
frontier bands (green) and the rest of the states (black). Up¬ 
per inset shows the first Brillouin zone with high symmetry 
points. 


to obtain tight-binding parameters [23-26]. This method 
has several advantages over the typical procedure of fit¬ 
ting a predetermined tight-binding model to the DFT 
bands; namely that the parameters relate to real elec¬ 
tronic distributions and values of the overlaps are inde¬ 
pendent of which hopping integrals one includes in the 
low-energy model. This approach to organic and organo- 
metallic molecular crystals, and the advantages this tech¬ 
nique has in these systems, is discussed in more detail 
elsewhere [23-25]. Most notably, molecular crystals tend 
to have well-separated bands, so the resulting WO basis 
is free from the difficulties of disentangling [26]. 

Fig. 2 shows the DFT band structure of Mo 3 S 7 (dmit) 3 . 
There are six isolated bands around the Fermi energy, all 
of which have hybrid Mo-dmit character. In the k x — k y 
plane two pairs of bands meet at the K and K' points (cf. 
Fig. 2, top left inset). Moving along k z , this degeneracy 
continues along the K — H and K' — H' lines. Near these 
lines the dispersion of these bands is linear in-plane, i.e., 
the in-plane mass vanishes. In the plane the dispersion is 
highly analogous to graphene, vida infra , and out of the 
plane the dispersion is approximately cosine-like. It is 
straightforward to show that, as for the honeycomb and 
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kagome lattices [27], the combination of time-reversal, in¬ 
version, and C 3 symmetries pin the Dirac points to the 
K — H lines. Only perturbations that break these sym¬ 
metries can move the Dirac point away from the K — H 
line or open a gap. 

The two bands that are not involved in the Dirac lines 
are extremely ID; with nearly-flat dispersion relations in 
the k z =0 (T — M — K) plane and cosine-like dispersion 
in the k z direction. The bands have slightly larger disper¬ 
sion in the k z = tt (H — L — A) plane, but even here the 
dispersion is very small compared to the in-plane hopping 
integrals (discussed below). 

Fig. 1(a) shows one of the WOs for Mo 3 S 7 (dmit) 3 . 
The six orbitals that contribute to the six low-energy 
bands are all related to this one by the C 3 symmetry of 
the molecules and inversion symmetry between the two 
molecules in a unit cell. These six orbitals are the frontier 
orbitals of the system, and are the basis of our models of 
its low energy physics. 

Note that along the T — A line the two lower bands 
are singly-degenerate, and the remaining bands form 
two doubly degenerate pairs. This is a straightforward 
consequence of the C 3 symmetry of the Mo 3 S 7 (dmit )3 
molecule. We label creation operators for the three sites 
of the triangle (WOs on the molecule) c\,i € {1,2,3} 
(henceforth we suppress spin labels). The molecular or¬ 
bital that belongs to the trivial (.A) irreducible repre¬ 
sentation is I'Fa) = ^(c{ + c\ + Cg)|0). The molecular 
orbitals that belong to the 2D E irreducible represen¬ 
tation may be written as I'I'b-) = ^(c{ — c|)| 0 ) and 

|4'e + ) = =(—c{ — c\ + 263 ) 10 ). We stress that, away 

from the T — A line, these labels are not good quantum 
numbers. Nevertheless, we will see below that this basis 
provides a simple physical description of the band struc¬ 
ture. 

To produce a low-energy model from the localized WO 
basis, one must truncate the set of all possible hopping 
integrals. The simplest method is to include all hopping 
integrals up to a real space cutoff, r c . Fig. 2 compares the 
DFT band structure with tight-binding models for r c = 
20 A and 35 A; the latter is essentially indistinguishable 
from the DFT bands. All parameters for these models 
are reported in the Supplementary Information [28]. It is 
worth noting that this model includes significant terms 
neglected by previous models of Mo 3 S 7 (dmit) 3 . 

Including only the largest three hopping integrals; 
tk = 59.69 eV, t g = 47.11 eV, and t z = 40.85 eV (twice 
as large as any other hopping integral); produces the 
kagomite model. This model captures all of the main 
features of the full DFT band structure (compare Figs. 
2 and 3), including the quasi-ID bands, and lets us ex¬ 
plain them. This model displays only cosine dispersion in 
the k z direction (see Fig. 3). Thus, all of the non-trivial 
phenomenology is confined to the plane and below we 
focus on the 2D kagomene lattice produced by setting 



FIG. 3: Main panel: Band structure of the kagomite lattice 
with tk, t g , and t z determined from the parameterization of 
Mo 3 S 7 (dmit )3 (cf. Table SI [28]). Note that the Dirac lines 
remain and the two ‘flat’ bands are now strictly one dimen¬ 
sional. The insets show the kagomene band structure for, 
from left to right, tk —» 00 , the DFT values, t g = 3tk/2, and 
t g —> 00 (in each case, the other t takes the DFT values and 
the filling is fixed at 2/3). It is clear from the t g —¥ 00 limit 
that the flat bands on the kagomene lattice are adiabatically 
connected to the flats bands in the two copies of the kagome 
lattice [29]. 


t z = 0. The origin of the kagomene lattice structure 
(Fig. 1(b)) of Mo 3 S 7 (dmit )3 can be straightforwardly 
understood in terms of the crystal structure and WO in 
Fig. 1(a). tk is the hopping between WOs on the same 
molecule whereas t g represents hopping between orbitals 
on adjacent molecules. 

The kagomene lattice interpolates between the honey¬ 
comb and kagome lattices. In the limit of t g —> ±00 
WOs on adjacent molecules are strongly hybridized and 
the model reduces to two decoupled kagome models in 
the basis of the bonding and anti-bonding hybrid or¬ 
bitals. Similarly in the molecular (tk —> ± 00 ) limit, the 
model is reduced to a honeycomb lattice formed from the 
A molecular orbitals, another honeycomb lattice formed 
from the E+ molecular orbitals, and two flat bands con¬ 
structed from the E — orbitals. The transition between 
the regimes exemplified by these two limits occurs at 
t g = 3tk/2 where there is an accidental three-fold de¬ 
generacy at F (see Fig. 3). Here the lower flat band 
switches which of the graphene-like bands it is degener¬ 
ate with. The DFT bands are in the molecular regime 
(t g < 3tk/2), consistent with intuition. 

To understand the flat bands in the kagomene band 
structure we introduce the ‘plaquette states’ (sketched 
in Fig. 1(b)); ^ E^i cos ( 777 ) 0 ] and = 

^ Sj =1 cos ( 7 rj /2 + 7 r/ 4 )cj, where j labels sequential 
sites around the dodecagon centered at R, such that sites 
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1 and 2 are on the same triangle. Explicit calculations 
shows that both of these states are eigenstates and are 
localised to a single dodecagon - this can be understood 
as follows. In both plaquette wavefunctions, amplitudes 
for the two vertices on the same dodecagon and the same 
triangle have the same magnitude and opposite sign, thus 
it is clear that they are, up to cyclic permutations, E— 
orbitals. The hopping terms tkc\c i and tkc\c 2 interfere 
destructively and cancel exactly. This applies to every 
triangle around the dodecagon. The two possible wave- 
functions on each ring are either bonding ( B or an¬ 
tibonding (AJ^) between adjacent triangles (Fig. 1(b)); 
hence they are split by an energy of 2 t g (cf. Fig. 3). 

These single dodecagon plaquettes can be combined to 
form extended plaquettes by taking linear combinations 
of contiguous dodecagon states. For open boundary con¬ 
ditions a sum over all dodecagons produces a plaquette 
state along the boundary of the system. However, for 
periodic boundary conditions, this sum vanishes since 
there is no boundary. Thus for a periodic system of 
N unit cells, which contains N dodecagons, there are 
2(N— 1) independent plaquette states. However, one can 
also construct four more localized states, which are de¬ 
generate with the plaquette states, by constructing topo¬ 
logically non-trivial superpositions of E— orbitals around 
the torus, Fig. 1(b). There are two distinct non-trivial 
loops and each may be either bonding or antibonding 
between triangles. Because of the non-trivial topology 
of these states, linear combinations with plaquette states 
can deform the loops but never close them. Now we have 

2 (TV + 1) degenerate states whereas each flat band con¬ 
tains N states. These extra states mean that there must 
be a point at which each of the flat bands touch one of 
the dispersive bands; this must occur at T [29]. 

Returning to kagomite we see that both the plaquette 
and topologically non-trival states move with a simple 
cosine dispersion in the z-direction. Additional inter¬ 
actions in the larger r c models mean that these states 
pick up a small in-plane dispersion, but it is clear that 
kagomite provides the basic framework to understand the 
DFT band structure. 

Finally, we turn to the role of interactions in this 
system. We have not attempted to calculate them di¬ 
rectly from our DFT calculations as this is known to 
be extremely difficult for tightly-packed molecular crys¬ 
tals such as these [23, 30-32]. However, it has recently 
been shown for Hubbard models on lattices formed of 
linked triangular molecules, such as the kagonrene and 
kagomite lattices, with n = 4 electrons per triangle that 
in the molecular ( t g —» 0) limit electronic correlations 
lead to the formation of localized spin-one moments with 
Heisenberg interactions [13, 14], In ID models this insu¬ 
lating phase has been shown to survive even far from the 
molecular limit [13, 14]. 

Each molecule is connected by three t z hopping inte¬ 
grals to the molecule above/below it, but only one t g 


hopping integral connects neighbours in the plane. It 
is straightforward, but tedious, to show that this means 
that as U —> oo, J z /J g —> 9 (t z /t g ) 2 . For Mo 3 S 7 (dmit) 3 , 
t z /tg = 0.468, so J z ~ 4.2 J g . This suggests that a ID 
spin-one model may be appropriate and that Haldane 
physics may be realised [13, 14], despite the intrinsically 
3D crystal structure and hopping integrals. 

If interchain interactions are relevant this would lead to 
a spin-one Heisenberg model on the layered honeycomb 
lattice. We are not aware of any studies of this model; 
however, as this lattice is bipartite one might expect that 
it will display long range order at sufficiently low tem¬ 
peratures if only nearest neighbour exchange interactions 
are included. However, there has recently been consid¬ 
erable interest in the possibility of spin liquid phases in 
the spin-1/2 J\ — J 2 model on the honeycomb lattice [33- 
35]. One is therefore tempted to speculate that second 
nearest neighbor exchange interactions could lead to a 
spin liquid state in MosS 7 (dmit) 3 . Interestingly, the ap¬ 
propriate next-nearest neighbour coupling introduces a 
chriality into the tight-binding model of Mo 3 S 7 (dmit) 3 , 
which is not unexpected for the P3 space group. It is 
important to note that while field theories of spin liq¬ 
uids are well known, it has proven harder to establish 
that realistic microscopic Hamiltonians in two or more 
dimensions have spin liquid ground states [7, 8, 33-38]. 
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